#ifndef _Cdd_External_Link_
#define _Cdd_External_Link_


#include <src/link/Qhull_Link.h>
#include <src/link/cdd/cdd_a.h>


namespace CDD {


	typedef particle<OS_double>										p_array_1;
	typedef particles_list<p_array_1>								p_array_2;
	typedef p_array_2												particles_list;
	typedef particles_list::Base									Base;
	typedef Base::value_type										particle_type;
	typedef Base::iterator											comp_iterator;
	typedef particle_type::iterator									part_iterator;
	typedef particles_list::indice_type								indice_type;
	typedef particles_list::indice_type								indice_type;
	typedef graph_base<particles_list>								graph_base_type;
	typedef graph_base_type::weight_set_type						weight_set_type;
	typedef graph_base_type::graph_type								edges_set_type;
	typedef graph_base_type::edge_type								edge_type;

	template<class A, class R>
	class cdd_edges_set_type : public QHULL::qh_edges_set_type<A,R>
	{
		template <class T>
		void construct()
		{
			OS_STATIC_CHECK(false); // you must implement the Delaunay graph for your own weight set type
		};
		template <>
		void construct<weight_set_default_def<base_value_type>>()
		{

			dd_set_global_constants();

			data_.clear();

			particle_type average = jmmath::moment1(*get_state_variable());
			dd_MatrixPtr M=NULL,M2=NULL,M3=NULL;
			dd_rowrange m_input = get_state_variable()->size()+1;
			dd_rowrange	i;
			dd_colrange d_input = get_state_variable()->dim()+2;
			dd_colrange j;
			dd_RepresentationType rep=dd_Generator;
			mytype value;
			OS_double rvalue;
			dd_rowset redrows,linrows,ignoredrows, basisrows;
			dd_colset ignoredcols, basiscols;

			dd_boolean found=dd_TRUE, newformat=dd_TRUE, successful=dd_TRUE, linearity=dd_FALSE;
			dd_NumberType NT = dd_Real;
			dd_ErrorType Error = dd_NoError;
			dd_SetFamilyPtr A=NULL;


			if (NT==dd_Unknown) {
			  Error = dd_ImproperInputFormat;
			  goto End_;
			}
			M=dd_CreateMatrix(m_input, d_input);
			M->representation=rep;
			M->numbtype=NT;
			OS_double renorm = 0;
			dd_set_d(value, 1.);
			for (i = 0; i < m_input; i++) {
				dd_set(M->matrix[i][d_input-1],value);
			}
			for (i = 0; i < m_input-1; i++) {
				particle_type local_particle((*get_state_variable())[i]-average);
				OS_double norm_2 = jmmath::norm2(local_particle);
				dd_set_d(value, norm_2);
				renorm = std::max(renorm,norm_2);
				dd_set(M->matrix[i][0],value);
				for (j = 1; j < d_input-1; j++) {
					dd_set_d(value, -2.*local_particle[j-1]);
					dd_set(M->matrix[i][j],value);
				}
			}
			dd_set_d(value, 1e+6);
			dd_set(M->matrix[m_input-1][0],value);
			dd_set_d(value, 0.);
			for (j = 1; j < d_input-1; j++) {
				dd_set(M->matrix[m_input-1][j],value);
			}
			
			A=dd_Matrix2Adjacency(M, &Error);
			if (Error != dd_NoError) goto End_;

			for (i=0; i<A->famsize; i++) {
				set_type  set = A->set[i];
				long elem;
				for (elem=i+1;elem<=set[0];elem++)
				{
					if (set_member(elem,set) && (elem!= m_input) ) {
						data_(i,elem-1)=1.;
						data_(elem-1,i)=1.;
					}
				}
			}
		End_ :
			dd_clear(value);

			if (A!=NULL) dd_FreeSetFamily(A);
			dd_FreeMatrix(M);
			if (Error!=dd_NoError) {
				dd_WriteErrorMessages(stderr,Error);
			}
		}
	};
};
#endif
